setwd(dirname(rstudioapi::getActiveDocumentContext()$path)) # set wd to source file location. For RStudio only
rm(list=ls())

library(tidyverse)
library(openxlsx)
library(lubridate)
library(zoo)
library(ggplot2)
library(dplyr)
library(texreg) 
library(MASS) 
library(dotwhisker) 
library(margins)
library(clubSandwich) 
library(sandwich)
library(lmtest) 
library(plm) 
library(pastecs) 
library(xtable)
library(pscl) 

options(max.print = 999999)
options(warn=-1)

source('functions.R')

data<-read.xlsx('analysis_dataset.xlsx',1) 
inflation_mentions<-read.xlsx('inflation_mentions.xlsx',1) 

### Items in dataset ###

# volume: Number of news articles on the BoE (month-newspaper)
# vol_lag: volume lagged by one month
# volume_ma: 5-point moving average of volume
# abs_dis_targ: Absolute distance of the relevant inflation rate from the relevant inflation target	
# abs_dis_lag1: abs_dis_targ lagged by one month
# cpi_yoy: Inflation rate (Consumer Prices Index (CPI) 12-month inflation rate (%))
# cpi_yoy_lag1: cpi_yoy lagged by one month
# letter: Dummy, value 1 when inflation deviated more than 1% from target and BoE governor had to sent letter to chancellor
# letter_lag1: letter lagged by one month
# unemploy: Monthly unemployment rate	
# unemp_lag1: unemploy lagged by one month
# gbp_usd: pound-dollar exchange rate (monthly average)
# gbp_usd_ch: change in gbp_usd (*100)
# gbpusd_log: absolute gbp_usd_ch, log transformed
# gdp_mom: month on month GDP growth	
# gdp_mom1: gdp_mom lagged by one month
# rate_lev: interest rate level	
# rate_absch: absolute change in rate_lev
# rate_anych: Dummy, value 1 when any change in interest rate took place
# bailouts: Dummy, value 1 when UK banks were bailed out
# qe_annc: Dummy, value 1 when QE programmes (QE1-5) was announced
# n_appnt: Number of new appointments announced (MPC, FPC, PRC)
# n_reappnt: Number of re-appointments announced (MPC, FPC, PRC)
# hearings: Number of parliamentary hearings of BoE officers	
# news: Number of news articles published by the BoE
# speeches: Number of speeches given by BoE officials	
# infl_rep: Dummy, value 1 when (quarterly) inflation report was published	
# brexit: Dummy, value 1 when Brexit-related events took place (see Appendix Table B4); relevant and previous month
# election: Dummy, value 1 when general election took place, relevant and previous month
# ukfsi: Level of systemic financial stress (UK Financial Stress Index (UKFSI))
# n_allappnt: Total number of appointment announcements (n_appnt + n_reappnt)
# log_allappnt: (1 + n_allappnt), log transformed  
# log_news: (1 + news), log transformed
# log_speeches: (1 + speeches), log transformed

# This creates a combined dataset where each unit is a month
data %>% group_by(month) %>%
  summarise(volume=sum(volume),
            vol_lag=sum(vol_lag),
            rate_absch=mean(rate_absch),
            qe_annc=mean(qe_annc),
            bailouts=mean(bailouts),
            abs_dis_lag1=mean(abs_dis_lag1),
            gbpusd_log=mean(gbpusd_log),
            unemp_lag1=mean(unemp_lag1),
            ukfsi=mean(ukfsi),
            hearings=mean(hearings),
            log_speeches=mean(log_speeches),
            log_allappnt=mean(log_allappnt),
            election=mean(election),
            brexit=mean(brexit),
            year=mean(year),
            infl_rep=mean(infl_rep))->combined

Sys.setlocale("LC_TIME", "C") # Language issue with zoo package for dates
data$month<-as.yearmon(data$month)
inflation_mentions$month<-as.yearmon(inflation_mentions$month)

cbPalette <- c('#000000', '#E69F00', '#56B4E9', '#009E73', '#F0E442', '#0072B2', '#D55E00', '#CC79A7')

### Figure 2 ### 

print(volume_fig<-ggplot(data,aes(x=month,y=volume_ma,group=newspaper,color=newspaper))+
  scale_colour_manual(values=cbPalette)+geom_line()+theme_classic()+ylab('')+labs(color='')+ggtitle('')+
  theme(legend.position='top',plot.title=element_text(hjust = 0.5,size=15),axis.text=element_text(size=15),
        legend.text = element_text(size = 15))+theme(axis.text.x=element_text(angle=45,hjust=1,colour='black'))+
  theme(axis.text.y = element_text(colour='black'))+scale_x_yearmon('',breaks=seq(1998,2020,2)))

### Table 1 - Main models ### 

# A. Policy decisions
mvol1z <- zeroinfl(volume ~ vol_lag + rate_absch + qe_annc + bailouts +
                     factor(infl_rep) + factor(newspaper) | factor(newspaper),
                   data=data, dist='negbin',link='logit')
mvol1zcl<-texreg::extract(clusteredSE(mvol1z,data,'newspaper'))

# B. Policy outcomes
mvol2z <- zeroinfl(volume ~ vol_lag + abs_dis_lag1 + gbpusd_log + unemp_lag1 + ukfsi + 
                     factor(infl_rep) + factor(newspaper) | factor(newspaper),
                   data=data, dist='negbin',link='logit')
mvol2zcl<-texreg::extract(clusteredSE(mvol2z,data,'newspaper'))

# C. Reinforced accountability
mvol3z <- zeroinfl(volume ~ vol_lag + hearings + log_speeches +
                     factor(infl_rep) + factor(newspaper) | factor(newspaper), # but why not appointments here?
                   data=data, dist='negbin',link='logit')
mvol3zcl<-texreg::extract(clusteredSE(mvol3z,data,'newspaper'))

# Full model, included controls
mvol4z <- zeroinfl(volume ~ vol_lag + rate_absch + qe_annc + bailouts +
                    abs_dis_lag1 + gbpusd_log + unemp_lag1 + ukfsi + 
                    hearings + log_speeches +
                    log_allappnt + election + brexit +
                    factor(infl_rep) + factor(newspaper) | factor(newspaper),
                  data=data, dist='negbin',link='logit')

mvol4zcl<-texreg::extract(clusteredSE(mvol4z,data,'newspaper'))

screenreg(list(mvol1zcl,mvol2zcl,mvol3zcl,mvol4zcl),stars=c(0.01,0.05,0.1))

### Figure 3 - Bank of England coverage, inflation mentions and distance from the target ### 

print(vol_inf_dist_fig<-ggplot(inflation_mentions,aes(x=month,y=volume_5ma,color="#56B4E9"))+
  geom_line(aes(y=volume_5ma,color="Number of articles"))+ 
  geom_line(aes(y=mentions_5ma,color="Occurrences of term \"inflation\""))+ 
  geom_line(aes(y=dist_target/.03,color="Distance from inflation target"))+ 
  theme_classic()+
  labs(color='')+
  ggtitle('')+
  scale_x_yearmon('',breaks=seq(1998,2020,2))+
  scale_y_continuous(name="Coverage of the Bank of England (5-point moving average)",
                     sec.axis = sec_axis(~.*.03, name="Distance from inflation target"))+
  scale_color_manual(values=c("Number of articles"="#56B4E9",
                              "Occurrences of term \"inflation\"" = "#E69F00",
                              "Distance from inflation target" = "#000000"))+
  theme(axis.text.x=element_text(angle=45,hjust=1,colour='black'),
        axis.text.y = element_text(colour='black'),
        legend.position='top',plot.title=element_text(hjust = 0.5,size=10),
        axis.text=element_text(size=10),
        legend.text = element_text(size = 10)))

### Figure 4 - Coefficient plot of Model 4 ### 
  
print(coefplot<-coef_zinb(mvol4zcl,drop=c('brexit','election','vol_lag','log_allappnt')) %>%
  by_2sd(data) %>%
  dwplot(ci=.9) %>%
  relabel_predictors(c(rate_absch='Abs. change bank rate',
                       qe_annc='QE announcement',
                       bailouts='Bank bailout',
                       abs_dis_lag1='Distance inflation target (t-1)',
                       gbpusd_log="Abs. change exchange rate (logged)",
                       unemp_lag1='Unemployment rate (t-1)',
                       ukfsi='Financial stress',
                       hearings='No. of hearings',
                       log_speeches='No. of speeches (logged)')) + 
  geom_vline(xintercept=0, colour="grey60",linetype=2) + xlab("Beta coefficient") + ylab("") + scale_x_continuous(limits = c(-.2,1)) +
  theme_bw()+theme(plot.title = element_text(size=15),axis.text=element_text(size=14),legend.title = element_blank(),legend.text = element_text(size = 14),
                   legend.position = 'off') + labs(title='DV: Coverage volume') +
  scale_color_manual(values = c('#000000','#56B4E9'),guide = guide_legend(ncol=2,reverse=T)))

################################################################################
################################### APPENDIX ###################################
################################################################################

####### DESCRIPTIVE STATISTICS - APPENDIX A ####### 

### Table A1 - Descriptive statistics ###

round(stat.desc(data$volume), 2) # Coverage volume
round(stat.desc(data$rate_absch), 2) # Rate change
round(stat.desc(data$qe_annc), 2) # QE announcement
round(stat.desc(data$bailouts), 2) # Bank bailout
round(stat.desc(data$abs_dis_targ), 2) # Distance from target
round(stat.desc(data$gbpusd_absch), 2) # Exchange rate change (original)
round(stat.desc(data$gbpusd_log), 2) # Exchange rate logged
round(stat.desc(data$unemploy), 2) # Unemployment rate
round(stat.desc(data$ukfsi), 2) # Financial stress
round(stat.desc(data$hearings), 2) # Hearings
round(stat.desc(data$speeches), 2) # Speeches (original) 
round(stat.desc(data$log_speeches), 2) # Speeches (log-transformed)
round(stat.desc(data$n_allappnt), 2) # Appointments (original) (sum new + reappoint)
round(stat.desc(data$log_allappnt), 2) # Appointments (log-transformed)
round(stat.desc(data$election), 2) # Elections
round(stat.desc(data$brexit), 2) # Brexit

### Table A2 - Correlation matrix ###

# Subsetting the relevant variables (in our main model)
data_corr <- data[, c('rate_absch','qe_annc','bailouts','abs_dis_lag1','gbpusd_log','unemp_lag1','ukfsi','hearings','log_speeches','log_allappnt','election','brexit')]
corstarsl(data_corr) # Used Pearson for all (at least for now)

corrmat <- xtable(corstarsl(data_corr), 
                  caption="Correlation matrix",
                  label="tab:corrmat", 
                  latex.environments="center")
print(corrmat)

### Table A3 - Correlation of coverage volume by newspaper ###

vol_bynp<-as.data.frame(data[data$newspaper=='The Times',]$volume)
vol_bynp$dm<-as.data.frame(data[data$newspaper=='Daily Mail',]$volume)
vol_bynp$grdn<-as.data.frame(data[data$newspaper=='The Guardian',]$volume)
colnames(vol_bynp)<-c('vol_times','vol_dm','vol_grdn')

corrmat_vol_bynp <- xtable(corstarsl(vol_bynp),caption="Correlation matrix",label="tab:corrmat_vol_bynp",latex.environments="center")
corstarsl(vol_bynp)

### Tables A5-A7 - Stationarity tests ###

pdata<-data.frame(split(data$volume,data$newspaper))

# Im-Pesaran-Shin test (2003) - Ha: at least one series is stationary 
# (errors can be serially correlated and serial corr can vary across cross sections; relies on ADF on each series)
purtest(pdata,lags='AIC',exo='intercept',test='ips') # reject null
purtest(pdata,lags='AIC',exo='trend',test='ips') # reject null

# Maddala-Wu test (1999) - Ha: at least one series is stationary
# (allows for different autoregressive coefficients among cross sections; relies on ADF on each series)
purtest(pdata,lags='AIC',exo='intercept',test='madwu') # reject null
purtest(pdata,lags='AIC',exo='trend',test='madwu') # reject null

# Levin-Lin-Chu test (2003) - Ha: all series are stationary 
# (assumption of homogeneous cross section, all have same autoregressive coefficient; relies on ADF on each series)
purtest(pdata,lags='AIC',exo='intercept',test='levinlin') # fails to reject null
purtest(pdata,lags='AIC',exo='trend',test='levinlin') # fails to reject null

### Figure A1 - Bank of England coverage vs inflation levels

print(vol_inf_lev_fig<-ggplot(inflation_mentions, aes(x=month,y=volume_5ma,color="#56B4E9")) + geom_line()+
        geom_line(aes(y=volume_5ma,color='Number of articles')) + 
        geom_line(aes(y=mentions_5ma,color='Occurrences of term \"inflation\"')) + 
        geom_line(aes(y=inflation/.05,color='Inflation rate')) + 
        theme_classic()+labs(color='')+ggtitle('')+
        theme(legend.position='top',plot.title=element_text(hjust = 0.5,size=10),axis.text=element_text(size=10),
              legend.text = element_text(size = 10))+theme(axis.text.x=element_text(angle=45,hjust=1,colour='black'))+
        theme(axis.text.y = element_text(colour='black'))+scale_x_yearmon('',breaks=seq(1998,2020,2))+
        scale_color_manual(values = c("Number of articles" = "#56B4E9","Occurrences of term \"inflation\"" = "#E69F00","Inflation rate" = "#000000"))+
        scale_y_continuous(name = "Coverage of the Bank of England (5-point moving average)",sec.axis = sec_axis(~.*.05, name="Inflation rate")))

####### ALTERNATIVE OPERATIONALISATIONS - APPENDIX B ####### 

### Table B1 - Any change in bank rate dummy instead of absolute amount of change ###

mvol1ratez<-zeroinfl(volume ~ vol_lag + rate_anych + qe_annc + bailouts + 
                       factor(infl_rep) + factor(newspaper) | factor(newspaper),
                     data=data, dist='negbin',link='logit')
mvol1ratezcl<-texreg::extract(clusteredSE(mvol1ratez,data,'newspaper'))

mvol4ratez<-zeroinfl(volume ~ vol_lag + rate_anych + qe_annc + bailouts + 
                       abs_dis_lag1 + gbpusd_log + unemp_lag1 + ukfsi + 
                       hearings + log_speeches +
                       log_allappnt + election + brexit +
                       factor(infl_rep) + factor(newspaper) | factor(newspaper),
                     data=data, dist='negbin',link='logit')
mvol4ratezcl<-texreg::extract(clusteredSE(mvol4ratez,data,'newspaper'))

screenreg(list(mvol1ratezcl,mvol2zcl,mvol3zcl,mvol4ratezcl),stars=c(0.01,0.05,0.1))

### Table B2 - Letter dummy instead of distance from inflation target

mvol2letz<-zeroinfl(volume ~ vol_lag + letter_lag1 + gbpusd_log + unemp_lag1 + ukfsi + 
                      factor(infl_rep) + factor(newspaper) | factor(newspaper),
                    data=data, dist='negbin',link='logit')
mvol2letzcl<-texreg::extract(clusteredSE(mvol2letz,data,'newspaper'))

mvol4letz<-zeroinfl(volume ~ vol_lag + rate_absch + qe_annc + bailouts + 
                      letter_lag1 + gbpusd_log + unemp_lag1 + ukfsi + 
                      hearings + log_speeches +
                      log_allappnt + election + brexit +
                      factor(infl_rep) + factor(newspaper) | factor(newspaper),
                    data=data, dist='negbin',link='logit')
mvol4letzcl<-texreg::extract(clusteredSE(mvol4letz,data,'newspaper'))

screenreg(list(mvol1zcl,mvol2letzcl,mvol3zcl,mvol4letzcl),stars=c(0.01,0.05,0.1))

### Table B3 - Inflation levels instead of distance from the target

mvol2inflz<-zeroinfl(volume ~ vol_lag + cpi_yoy_lag1 + gbpusd_log + unemp_lag1 + ukfsi + 
                       factor(infl_rep) + factor(newspaper) | factor(newspaper),
                     data=data, dist='negbin',link='logit')
mvol2inflzcl<-texreg::extract(clusteredSE(mvol2inflz,data,'newspaper'))


mvol4inflz<-zeroinfl(volume ~ vol_lag + rate_absch + qe_annc + bailouts + 
                       cpi_yoy_lag1 + gbpusd_log + unemp_lag1 + ukfsi + 
                       hearings + log_speeches +
                       log_allappnt + election + brexit +
                       factor(infl_rep) + factor(newspaper) | factor(newspaper),
                     data=data, dist='negbin',link='logit')
mvol4inflzcl<-texreg::extract(clusteredSE(mvol4inflz,data,'newspaper'))

screenreg(list(mvol1zcl,mvol2inflzcl,mvol3zcl,mvol4inflzcl),stars=c(0.01,0.05,0.1))

### Table B4 - GBP/USD exchange rate instead of absolute change

mvol2gbpz <- zeroinfl(volume ~ vol_lag + abs_dis_lag1 + gbp_usd_ch + unemp_lag1 + ukfsi + 
                     factor(infl_rep) + factor(newspaper) | factor(newspaper),
                   data=data, dist='negbin',link='logit')
mvol2gbpzcl<-texreg::extract(clusteredSE(mvol2gbpz,data,'newspaper'))

mvol4gbpz <- zeroinfl(volume ~ vol_lag + rate_absch + qe_annc + bailouts +
                     abs_dis_lag1 + gbp_usd_ch + unemp_lag1 + ukfsi + 
                     hearings + log_speeches +
                     log_allappnt + election + brexit +
                     factor(infl_rep) + factor(newspaper) | factor(newspaper),
                   data=data, dist='negbin',link='logit')

mvol4gbpzcl<-texreg::extract(clusteredSE(mvol4gbpz,data,'newspaper'))

screenreg(list(mvol1zcl,mvol2gbpzcl,mvol3zcl,mvol4gbpzcl),stars=c(0.01,0.05,0.1))

### Table B5 - Month-on-month GDP growth instead of unemployment rate ###

mvol2gdpz<-zeroinfl(volume ~ vol_lag + abs_dis_lag1 + gbpusd_log + gdp_mom1 + ukfsi + 
                factor(infl_rep) + factor(newspaper) | factor(newspaper),
                data=data, dist='negbin',link='logit')
mvol2gdpzcl<-texreg::extract(clusteredSE(mvol2gdpz,data,'newspaper'))


mvol4gdpz<-zeroinfl(volume ~ vol_lag + rate_absch + qe_annc + bailouts + 
                abs_dis_lag1 + gbpusd_log + gdp_mom1 + ukfsi + 
                hearings + log_speeches +
                log_allappnt + election + brexit +
                factor(infl_rep) + factor(newspaper) | factor(newspaper),
                data=data, dist='negbin',link='logit')
mvol4gdpzcl<-texreg::extract(clusteredSE(mvol4gdpz,data,'newspaper'))

screenreg(list(mvol1zcl,mvol2gdpzcl,mvol3zcl,mvol4gdpzcl),stars=c(0.01,0.05,0.1))

### Table B6 - UKFSI without bailout and QE dummy ###

mvol4ukfisz<-zeroinfl(volume ~ vol_lag + rate_absch + 
                       abs_dis_lag1 + gbpusd_log + unemp_lag1 + ukfsi + 
                       hearings + log_speeches +
                       log_allappnt + election + brexit +
                       factor(infl_rep) + factor(newspaper) | factor(newspaper),
                     data=data, dist='negbin',link='logit')
mvol4ukfiszcl<-texreg::extract(clusteredSE(mvol4ukfisz,data,'newspaper'))

screenreg(list(mvol1zcl,mvol2zcl,mvol3zcl,mvol4ukfiszcl),stars=c(0.01,0.05,0.1))

### Table B7 - News releases instead of speeches ###

mvol3newsz<-zeroinfl(volume ~ vol_lag + hearings + log_news +
                factor(infl_rep) + factor(newspaper) | factor(newspaper),
                data=data, dist='negbin',link='logit')
mvol3newszcl<-texreg::extract(clusteredSE(mvol3newsz,data,'newspaper'))

mvol4newsz<-zeroinfl(volume ~ vol_lag + rate_absch + qe_annc + bailouts + 
                abs_dis_lag1 + gbpusd_log + unemp_lag1 + ukfsi + 
                hearings + log_news +
                log_allappnt + election + brexit +
                factor(infl_rep) + factor(newspaper) | factor(newspaper),
                data=data, dist='negbin',link='logit')
mvol4newszcl<-texreg::extract(clusteredSE(mvol4newsz,data,'newspaper'))

screenreg(list(mvol1zcl,mvol2zcl,mvol3newszcl,mvol4newszcl),stars=c(0.01,0.05,0.1))

####### ALTERNATIVE SPECIFICATIONS - APPENDIX C ####### 

### Table C1 - Main models with linear time trend ###

mvol1ztime <- zeroinfl(volume ~ vol_lag + rate_absch + qe_annc + bailouts +
                         factor(infl_rep) + factor(newspaper) + as.integer(year-1997) | factor(newspaper),
                       data=data, dist='negbin',link='logit') 
mvol1zcltime<-texreg::extract(clusteredSE(mvol1ztime,data,'newspaper'))

mvol2ztime <- zeroinfl(volume ~ vol_lag + abs_dis_lag1 + gbpusd_log + unemp_lag1 + ukfsi + 
                         factor(infl_rep) + factor(newspaper) + as.integer(year-1997)  | factor(newspaper),
                       data=data, dist='negbin',link='logit')
mvol2zcltime<-texreg::extract(clusteredSE(mvol2ztime,data,'newspaper'))

mvol3ztime <- zeroinfl(volume ~ vol_lag + hearings + log_speeches +
                         factor(infl_rep) + factor(newspaper) + as.integer(year-1997)  | factor(newspaper), # but why not appointments here?
                       data=data, dist='negbin',link='logit')
mvol3zcltime<-texreg::extract(clusteredSE(mvol3ztime,data,'newspaper'))

mvol4ztime <- zeroinfl(volume ~ vol_lag + rate_absch + qe_annc + bailouts +
                         abs_dis_lag1 + gbpusd_log + unemp_lag1 + ukfsi + 
                         hearings + log_speeches +
                         log_allappnt + election + brexit +
                         factor(infl_rep) + factor(newspaper) + as.integer(year-1997)  | factor(newspaper),
                       data=data, dist='negbin',link='logit')

mvol4zcltime<-texreg::extract(clusteredSE(mvol4ztime,data,'newspaper'))

screenreg(list(mvol1zcltime,mvol2zcltime,mvol3zcltime,mvol4zcltime),stars=c(0.01,0.05,0.1))

### Table C2 - Main models with normalised dependent variable ###

data %>% group_by(newspaper) %>% mutate(vol_norm = (volume-min(volume))/(max(volume)-min(volume)),
                                        vol_lag_norm = (vol_lag-min(vol_lag,na.rm=T))/(max(vol_lag,na.rm=T)-min(vol_lag,na.rm=T)))->data

mvol1norm<-lm(vol_norm ~ vol_lag_norm + rate_absch + qe_annc + bailouts + 
                factor(infl_rep) + factor(newspaper),data=data)
mvol1normcl<-texreg::extract(coeftest(mvol1norm, vcov = vcovCL, cluster = ~newspaper))

mvol2norm<-lm(vol_norm ~ vol_lag_norm + abs_dis_lag1 + gbpusd_log + unemp_lag1 + ukfsi +
                factor(infl_rep) + factor(newspaper),data=data)
mvol2normcl<-texreg::extract(coeftest(mvol2norm, vcov = vcovCL, cluster = ~newspaper))

mvol3norm<-lm(vol_norm ~ vol_lag_norm + hearings + log_speeches +
                factor(infl_rep) + factor(newspaper),data=data)
mvol3normcl<-texreg::extract(coeftest(mvol3norm, vcov = vcovCL, cluster = ~newspaper))

mvol4norm<-lm(vol_norm ~ vol_lag_norm + rate_absch + qe_annc + bailouts + 
                abs_dis_lag1 + gbpusd_log + unemp_lag1 + ukfsi +
                hearings + log_speeches +
                log_allappnt + election + brexit +
                factor(infl_rep) + factor(newspaper),data=data)
mvol4normcl<-texreg::extract(coeftest(mvol4norm, vcov = vcovCL, cluster = ~newspaper))

screenreg(list(mvol1normcl,mvol2normcl,mvol3normcl,mvol4normcl),stars=c(0.01,0.05,0.1))

### Table C3 - Negative binomial models ### 

mvol1<-glm.nb(volume ~ vol_lag + rate_absch + qe_annc + bailouts + 
                factor(infl_rep) + factor(newspaper),data=data)
mvol1cl<-texreg::extract(coeftest(mvol1, vcov = vcovCL, cluster = ~newspaper))

mvol2<-glm.nb(volume ~ vol_lag + abs_dis_lag1 + gbpusd_log + unemp_lag1 + ukfsi +
                factor(infl_rep) + factor(newspaper),data=data)
mvol2cl<-texreg::extract(coeftest(mvol2, vcov = vcovCL, cluster = ~newspaper))

mvol3<-glm.nb(volume ~ vol_lag + hearings + log_speeches +
                factor(infl_rep) + factor(newspaper),data=data)
mvol3cl<-texreg::extract(coeftest(mvol3, vcov = vcovCL, cluster = ~newspaper))

mvol4<-glm.nb(volume ~ vol_lag + rate_absch + qe_annc + bailouts + 
                abs_dis_lag1 + gbpusd_log + unemp_lag1 + ukfsi +
                hearings + log_speeches +
                log_allappnt + election + brexit +
                factor(infl_rep) + factor(newspaper),data=data)
mvol4cl<-texreg::extract(coeftest(mvol4, vcov = vcovCL, cluster = ~newspaper))

screenreg(list(mvol1cl,mvol2cl,mvol3cl,mvol4cl),stars=c(0.01,0.05,0.1))

### Table C4 - Vuong test for appropriateness of zero-inflation ###

vuong(mvol4,mvol4z) # we can reject the null that models are indistiguishible: table indicates which model is better (in our case model 2, zeroinflated)

### Table C5 - OLS models ###

mvol1ols<-lm(volume ~ vol_lag + rate_absch + qe_annc + bailouts + 
               factor(infl_rep) + factor(newspaper),data=data)
mvol1olscl<-texreg::extract(coeftest(mvol1ols, vcov = vcovCL, cluster = ~newspaper))

mvol2ols<-lm(volume ~ vol_lag + abs_dis_lag1 + gbpusd_log + unemp_lag1 + ukfsi +
               factor(infl_rep) + factor(newspaper),data=data)
mvol2olscl<-texreg::extract(coeftest(mvol2ols, vcov = vcovCL, cluster = ~newspaper))

mvol3ols<-lm(volume ~ vol_lag + hearings + log_speeches +
               factor(infl_rep) + factor(newspaper),data=data)
mvol3olscl<-texreg::extract(coeftest(mvol3ols, vcov = vcovCL, cluster = ~newspaper))

mvol4ols<-lm(volume ~ vol_lag + rate_absch + qe_annc + bailouts + 
               abs_dis_lag1 + gbpusd_log + unemp_lag1 + ukfsi +
               hearings + log_speeches +
               log_allappnt + election + brexit +
               factor(infl_rep) + factor(newspaper),data=data)
mvol4olscl<-texreg::extract(coeftest(mvol4ols, vcov = vcovCL, cluster = ~newspaper))

### Table C6 - Main models with combined volume variable ###

mvol1comb<-glm.nb(volume ~ vol_lag + rate_absch + qe_annc + bailouts + 
                    factor(infl_rep),data=combined)

mvol2comb<-glm.nb(volume ~ vol_lag + abs_dis_lag1 + gbpusd_log + unemp_lag1 + ukfsi +
                    factor(infl_rep),data=combined)

mvol3comb<-glm.nb(volume ~ vol_lag + hearings + log_speeches +
                    factor(infl_rep),data=combined)

mvol4comb<-glm.nb(volume ~ vol_lag + rate_absch + qe_annc + bailouts + 
                    abs_dis_lag1 + gbpusd_log + unemp_lag1 + ukfsi +
                    hearings + log_speeches +
                    log_allappnt + election + brexit +
                    factor(infl_rep),data=combined)

screenreg(list(mvol1comb,mvol2comb,mvol3comb,mvol4comb),stars=c(0.01,0.05,0.1))

####### MODELS BY NEWSPAPER - APPENDIX D ####### 

### Table D1 - Models by newspaper (negative binomial) ###

mvol4times<-glm.nb(volume ~ vol_lag + rate_absch + qe_annc + bailouts + 
                abs_dis_lag1 + gbpusd_log + unemp_lag1 + ukfsi + 
                hearings + log_speeches +
                log_allappnt + election + brexit +
                factor(infl_rep),data=data[data$newspaper=='The Times',])

mvol4grdn<-glm.nb(volume ~ vol_lag + rate_absch + qe_annc + bailouts + 
                     abs_dis_lag1 + gbpusd_log + unemp_lag1 + ukfsi + 
                     hearings + log_speeches +
                     log_allappnt + election + brexit +
                     factor(infl_rep),data=data[data$newspaper=='The Guardian',])

mvol4mail<-glm.nb(volume ~ vol_lag + rate_absch + qe_annc + bailouts + 
                     abs_dis_lag1 + gbpusd_log + unemp_lag1 + ukfsi + 
                     hearings + log_speeches +
                     log_allappnt + election + brexit +
                     factor(infl_rep),data=data[data$newspaper=='Daily Mail',])

screenreg(list(mvol4times,mvol4grdn,mvol4mail),stars=c(0.01,0.05,0.1))

